SEMIGROUP SPLITTING AND CUBATURE APPROXIMATIONS 
FOR THE STOCHASTIC NAVIER-STOKES EQUATIONS* 

PHILIPP DORSEKt 

Abstract. Approximation of the marginal distribution of the solution of the stochastic Navier- 
Stokes equations on the two-dimensional torus by high order numerical methods is considered. The 
corresponding rates of convergence are obtained for a splitting scheme and the method of cubature on 
Wiener space applied to a spectral Galerkin discretisation of degree N. While the estimates exhibit 
a strong A'^ dependence, convergence is obtained for appropriately chosen time step sizes. Results of 
numerical simulations are provided, and confirm the applicability of the methods. 
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1. Introduction. The issue of turbulence in fluid flows is an essentially unsolved 
', problem. From the perspective of numerical analysis, its main difficulty is that a 

' direct numerical simulation (DNS), resolving all relevant temporal and spatial scales, 

, is unavailable for many practically relevant geometries. Hence, we can only use results 

' from underresolved simulations, which are often useless due to their severely reduced 

accuracy. 

This has led to the introduction of reduced models that deal with the closure 
problem, see e.g. |30l These models deal with underresolution by introducing a 
model for the effects taking place on scales smaller than those that are resolved. 

We are concerned with a different approach to turbulence modelling. In the last 
years, the introduction of noise into the equations of fluid dynamics has become the 
focus of research (see e.g. [H HU [SI [371 [T] ) . In particular, Hairer and Mattingly proved 
ly-^ ' in [13l [14] that the stochastic Navier-Stokes equations on the two-dimensional torus 

OA i with finite-dimensional, additive noise have ergodic dynamics, and estimated the rate 

' of convergence to the invariant measure. 

, There are several ways to find numerical approximations to stochastic differen- 

tial equations. Pathwise and strong schemes aim to find path-by-path simulations 
of the problem, aiming for convergence almost surely or in mean square error; see 
e.g. [m [19] for recent survey articles. In this work, we consider the approximation 
ij , of the stochastic Navier-Stokes equations on the two-dimensional torus by weak ap- 

I^S ' proximation schemes. As we are actually searching for estimates of functionals of 

, the invariant measure, this allows us to benefit from the typically higher rate of con- 

vergence of weak schemes. Related work concerning the numerical analysis of weak 
methods for stochastic partial differential equations can be found in [S] [TU] . 

In contrast to [16], we propose a simulation method, based either on a splitting 
similar to the Ninomiya-Victoir scheme [29], or the method of cubature on Wiener 
space of Kusuoka [22] and Lyons and Victoir [23] , extending the applicability of such 
approaches from bounded vector fields to vector fields that are neither Lipschitz con- 
tinuous nor linearly bounded anymore (see also [2] for other such extensions, in that 
case to vector fields behaving like square roots at zero) . The advantage of such an ap- 
proach is that it is trivial to parallelise, as every path can be simulated independently. 
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Furthermore, in the case of sphtting schemes, we can furthermore reuse weU-tested, 
robust and fast solvers for the deterministic Navier-Stokes or Euler equations to obtain 
solvers for the stochastic Navier-Stokes equations with minimal effort. 

To derive rates of convergence, we employ the theory derived in |10j . While we 
are unable to prove rates of convergence on the continuous level, a discretisation by 
a spectral Galerkin scheme allows us to obtain a suboptimal convergence estimate. 

This paper is organised as follows. In Section [21 we recall the definition of the 
stochastic Navier-Stokes equations in the setting of Hairer and Mattingly and consider 
them from the perspective of weighted spaces used in TU' . Our analysis profits greatly 
from the fundamental results shown by Hairer and Mattingly in j24j [131 [M] . Section[3] 
is devoted to the derivation of estimates for the error done by a spectral Galerkin 
approximation. Section [?] presents the main results of this paper, estimates for full 
discretisations of the stochastic Navier-Stokes equations by splitting and cubature 
schemes. In Section [51 we present the results of numerical calculations for a model 
problem with ergodic dynamics, and in Section [6l we sum up our results. 

2. The stochastic Navier-Stokes equations and vifeighted spaces. Con- 
sider, as in [13l[T4], the vorticity formulation of the stochastic Navier-Stokes equations 
on the two-dimensional torus T^, 

d 

dw{t,wo) ^ iyAw{t,wo)dt + B{ICw{t,wo),w{t,wo))dt + ^qjfk,dWj' , (2.1) 
w{0,wo) = Wo- 

The state space is L^, the space of mean zero square integrable functions, with norm 
ll'll and scalar product (•, •). Furthermore, A is the Laplacian, /C the inverse of the 
rotation VAu — d2Ui — diU2 in the space of divergence free vector fields, VA(/Cw) = w 
and V -JCw — 0, i?(u, w) — ~{u-V)w the Navier-Stokes nonlinearity, and {Wl)j=i^... ^j. 
a d-dimensional Brownian motion. The qj are nonvanishing real numbers, qj G M\{0}, 
and fk are the orthonormal eigenfunctions of A on T^, 

'~\{27r^)-^/^cos{k-x), else, ^ ' 

where 

Zl := {k ^ (fci, fcs) e : either fca > 0, or fca = and ki > 0} . (2.3) 

We also define the Sobolev spaces of divergence-free, mean zero functions with 



norm EfeGZ^^fc/felU '■= y Z]/cez2(fci + ^D'^l^fcP' which is non-degenerate due to 
the mean zero condition (the term for k — (0,0) vanishes). We note, in particular, 
that 

- {Aw,w) ^\\w\\l. (2.4) 

Similarly as in 14, Section 5.3], we introduce the weight iprjiw) :— exp(77||u'||^) 
with some 77 > and consider the weighted space B'^'^{L,'^), given as the closure of 
the space of smooth, cylindrical functions / — .g((-, ei), . . . , (•, e„)), g G C|f (R"), with 
respect to the norm 

WfU, ■■= sup ^P,iw)-'\f{w)\. (2.5) 
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For a more complete exposition of such spaces, we refer the reader to [10] . We remark 
that [ini Section 5] provides several conditions under which Markov semigroups of 
operators defined on such weighted spaces become strongly continuous, generalising 
the ideas in [Ml Section 5.3]. 

Proposition 2.1. The Markov semigroup {Pt)t>o defined through Ptf{wo) := 
K[f{w{t,wo))] is strongly continuous on yB^''(L^) for 77 > small enough. 

Proof. This follows from lOj Theorem 5.8] and [131 Theorem A. 3]. A very similar 
result is proved in [14, Theorem 5.10]. □ 

Contrary to the approach used in [29l |32] , we are not able to split this problem 
into a part corresponding fully to the drift and another for the diffusion: the pro- 
cess y{t,WQ)t :— wq + J2j=i Ij^t fkj corresponding to the diffusion does not satisfy 
E[7/'^(2/(i, Wo))] < Ktpriiwo) with K > constant for t small enough, which means 
that we cannot use standard Ninomiya-Victoir splittings. 

Thus, we split up the equation differently. For a given e e (0, 1), we introduce 
the deterministic vorticity equation, 

^w\t,wo) = {l-e)vl\w\t,wo) + B{lCw\t,wo),w^{t,wo)), (2.6) 
at 

w^(0, Wo) = Wo, 

and a stochastic heat equation defining an Ornstein-Uhlenbeck process on L'^ , 

d 

du;2(i,u;o) -£i^Aw2(t,u;o)di + ^gj/fc^.dVF/, w'^{0,wq) ^ wq. (2.7) 

Define by Plf{wo) := E[f{w^{t,wo))] and P^f{wo) := E[f{w^{t,wo))] the Markov 
semigroups corresponding to and uj^. 

Lemma 2.2. For 77 > 0, (P/)f>o defines a strongly continuous semigroup on 
B^-ih^) with ||P/IL(e..(L2)) < 1. 

Proof. The strong continuity is obtained using (10^ Theorem 5.8]. The necessary 
bounds are proved by applying [TH Theorem A. 3]; see also [131 Theorem 5.10]. 

The deterministic vorticity equations have L^-contractive dynamics, as 

\\w^{t,wo)\\'^ = ||wo||^ + / (eJ^Au;^(s, Wo) + B(/Cw^ (s, wq), w^(s, wo)), w^(s, wo))ds 
Jo 

< llw^olP, (2.8) 

which yields the norm bound. The proof is thus complete. □ 

The cumbersome proof of the following proposition is postponed to the appendix. 
Proposition 2.3. If r] > is small enough, there exists uj > such that the 

process 1 1— exp(— wt)7/'^(w^(i, wo)) is a positive supermartingale, i.e. 

E[V'r,(w2(i,wo)] < exp(a;i)V'r,(w'^(i,«^o)). (2.9) 

Lemma 2.4. For rj > small enough, {P^)t>Q is strongly continuous on B'^''' {iJ^) 
with bound \\Pt\\L{B''''^(L^)) — exp(wt). 

Proof. Clear from Proposition [53] (see also [IHl Example 5.4]). □ 
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3. Spectral Galerkin approximations for stochastic Navier-Stokes equa- 
tions. For the stochastic Navier-Stokes equations, we cannot argue directly as in [TU] : 
there do not appear to be useful weight functions on spaces of more regular functions 
(such spaces are nevertheless invariant with respect to the dynamics of (|2.ip : see [24l 
Section 3.4] in this regard). We will therefore settle with a weaker result: we shall 
prove that spectral Galerkin approximations using Fourier modes up to degree N 
yield a convergent scheme, which can then be approximated by a splitting or a cu- 
bature scheme with A^-dependent error bound. As the A^-dependence of the estimate 
is given explicitly, we can derive convergent schemes by choosing the time step size 
small enough in relation to N. 

Consider therefore the spectral Galerkin approximation of (|2.1I) . 

dwAr(t, Wo) = J^Awjv(i, wo)di (3.1a) 



+ TTNB{ICwN{t,wo),WN{t,wo))dt + y^gj/fc,dW^/, 

WAr(0,wo) = ttatWo, (3.1b) 

see also [TT], where ttn: — > is the projection onto the space Hn of tensor 
products of trigonometric polynomials of degree N, 

Un := spanj/fc: max|A:i| < 7v| , (3.2) 

and N is assumed to be large enough so that fk^ G Ti-N for j = 1, . . . Its split 
semigroups are given by 

^w]^{t,WQ) = TTNB{JCw\i{t,'W()),w\i{t,'W())), w]v(0, Wo) = Wo, and (3.3) 

d 

dw%{t,wo) ^ uAw%{t,WQ)dt + y^^qjfk^dW^ , w%{0,wq)=wq. (3.4) 

3 = i 

The choice s = 1 made here is not admissible above: in the space continuous setting, 
the results from |T4] do not allow us to apply [lOl Theorem 5.8] to conclude that 

is strongly continuous for this choice. As Hn is finite-dimensional, however, 
we do not have to distinguish between different topologies, and it follows that the 
Markov semigroups , P-t^'^ and P^'^ of wat, w\j and w^^ are strongly continuous 
on B^^{'Hn) if 77 > is small enough. In case that a solver for deterministic Navier- 
Stokes equations is available, it is also possible to use e < 1 here (the case e = 1 
corresponds to splitting up into a deterministic Euler equation). 

We now estimate the error of the spectral Galerkin approximation. 

Proposition 3.1. For any a > 0, Wq e and t > 0, 

||w(f,wo)-wjv(i,wo)|p < CN-^\\w{t,wo)\\l 

+ C^N-'exp(c^t+^ ||w(a,wo)||?da^^ \\w{s,wo)\\tds. (3.5) 

Proof. Let eAr(t) :— 7rArw(i, wo) — WAr(t, wo) G Hn and r]N{t) := w(t,wo) — 
TrN'w{t,wo). Then, 

deN{t) — i^AeN{t) + ttn {B{JCwN[t, wo), eN{t)) + B{ICeN{t),TTNw{t, wo))) dt 

+ TTN (B(JCTiNw(t, Wo), riN(t)) + B{ICr]N{t), w{t, wq))) dt. (3.6) 
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It results that 

^;^l|eivWII' = -i^\\eNit)\\ldt + {BilCeNit),TTNw{t,wo)),eNit)) 

+ {B{ICTrN'w{t,wo),riNit)) + B{ICT]N{t),w{t,wa)),eNit)). (3.7) 

We now proceed similarly as in [13, Proof of Lemma 4.10, point 3]. For any S > 0, 
we estimate 

\{B{ICh,w),0\ < S\m + ^llCir + ^Ml\\h\\'. (3.8) 

This yields 

\{B{ICeN{t),TrNwit,wo)),eNit))\ < <5||e^v(t)||? + ^lle^vWH' 



+ -\\TrNw{t,wa)\\l\\eN{t)f and (3.9) 



-\\w{t,wo)\\l\\mit)r. (3.10) 



For the final term, we apply 



m]ch,w)x)\<s\\c\\i + ^^\\h\\i\\wr, (3.11) 



which shows 



\{B{ICTTNwit,wo),VN{t)),eNit))\ < S\\eNit)\\l + Slkjvw(t, wo)li?|| W Wf • (3.12) 

Ad 



Choosing S — ^ and combining the above estimates yields 

^^||eAr(t)|p < -"^WeNiml + ^lleA'WII' + f II W(t, wo)||?||eAr(i)|" 
a,, / ^,,9 3C 



+ i^-\\w{t,wom + —\\TTNw{t,wo)\\ij llwWII • (3.13) 

Using IIttjv'il'IIi < ||w||i, we obtain 

l^^W^NitW < {Cc + ^\\w{t,wo)\\l) IWe^l' + C^\\w{t,wo)\\l\\rMtW- (3-14) 
An application of Gronwall's inequality yields, as eAr(O) = 0, 



(3.15) 



^\\eNit)r< I CJwis,wo)\\i\\vNis)ry< 



X exp ^Ca{t ^ ^) + ^ ^ W^i'^^ ^o)\\i<ia^ ds. 



As ||ti; — nffw\\ < CN we see that ||77Ar(t)|| < CN ^\\w{t,wo)\\i, which yields 

l\\eN{t)f < CaN-^ \\w{s,wo)\\tcxp (Ca.it - s) + ^ \\w{a,w„)\\lda^ ds 

<Cc.N-^exp(cc,t+^ f \\w{a,wo)\\lda\ f \\w{s,wo)\\tds. (3.16) 
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The result follows due to ||w(t,Wo) — WAr(t,u'o)|| < ||eAr(t)|| +CA^^-^||u)(t,u)o)||i. □ 
Corollary 3.2. For any wo € and T > 0, there exists a constant C = 
Cwo,T > such that for any t £ [0, T], 



¥.[\\w{t,WQ)-WN{t,WQ)f] <CN-\ 



(3.17) 



Proof. From ProDOsition l3.1l and an application of the Cauchy-Schwarz inequality, 
we see that we need to prove 



E[\\w{t,wo)\\l]+E 



+ E 



cxp ( a I ||w(cr, u;o)||idCT 

2 



\w{s,wo)\\ids 



< K 



(3.18) 



for all t G [0, T] with some K = Kt^wo > 0. For the first and third term, this follows 
from '24| Theorem 3.7], and for the second, from |13j Lemma 4.10]. □ 

Remark 3.3. Actually, it seems quite plausible here that the assumption wq G 
is too strong. Indeed, the results in i25f show that if Wq G L^, then w{t,WQ) € H** 
for all s > for subsequent times, and \2b\ Lemma A. 3] gives some quantitative 
estimates. It remains unclear to us however how this can be used to prove an estimate 



forE 



f^\\w{s,wo)\\idi 



The estimate from Corollary 13.21 allows us to estimate the pointwise approxima- 
tion error of the weak approximation of the stochastic Navier-Stokes equation by the 
spectral Galerkin scheme. 

Theorem 3.4. Assume ip e F^'"(L2) n C^{h'^) with 



Cip :— sup ipfj{w) ^\\Dlp{w)\\ < oo 



(3.19) 



for some fj G [0,r]/2]. Then, for w G and T > 0, there exists a constant C = Cw,T,ip 
such that for all t G [0, T], 



\PMy^)~P!'Mn.)H\<CN-\ 



(3.20) 



Proof. By the fundamental theorem of calculus, 

\(p{w{t, Wo)) - (p{wN{t, Wo))\ (3.21) 

< / \\Dp{ew{t,wo) + {l-e)wN{t,wo))\\-\\w{t,wo)-WN{t,wo))\\d0. 
Jo 

The assumption on (p together with the convexity of w h- !■ exp(77||?i'|P) yields 
\\Dip{0wit, wo) + (1 - 0)wNit, wo))\\ 

< (exp(?7||w(t,wo)|P) +cxp(fy||wAr(t,wo)|P)) • 
Therefore, the Cauchy-Schwarz inequality implies 

\PMw) ~ F,^(^|«„(u;)| < C^E [\\w{t,wo) - WN{t,wo)\\^]^^^ x 

X (E[cxp{2i^\\wit,wo)f)]'^^ +E[eM'2ij\\wNit,wo)f)]'^^] 



(3.22) 



(3.23) 
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Note that the estimate in [T31 Lemma 4.10, 1.] also holds true for wjv(t, Wq) instead 
oiw{t,wo)- Therefore. Corollary 13.21 proves the claimed estimate. □ 

In the discrete setting, it is easy to analyse the differential operators corresponding 
to the split semigroups. For A; > 0, we define B^'\'Hn) as the closure of C'^{'Hn) 
with respect to the norm 

k 

ll/IU„.fe:=ll/ll0.+El/k.^' (3.24) 
i=i 

where the seminorms j = 1, . . . , fc, are given by 

\f\i>vJ sup V'^(w)"^||£'^/(w)|L((«„)«..;E). (3.25) 

We denote by Gf with domain domQ^ the infinitesimal generator of (P/^'"')t>o, 
j — 1, 2, and by with domain domQ^ the infinitesimal generator of {Pf^)t>o. 
Lemma 3.5. For any e > 0, 

Bp {Hn) C dom n dom n dom . (3.26) 

For k > 0, , gf : Bfl^inN) Bp+'inw), j = 1,2, are continuous operators, 
and 

Furthermore, 

g^ip = g^ip + g^ip for all ip e Bp (n^). (3.28) 



Proof. For ip G B'j^^^i'^N), we see by the fundamental theorem of calculus and 
the estimates in [Mj Appendix] that with a > 0, 

|gf^(u;)| = \D^{w) iTTNBilCw,w))\ < \\Dp:iw)\\ ■ {N^+^w\f) 

< CiV^exp(£|lu;||2)||i:>(^(u))||, (3.29) 



1 

g^ifiiw)] = |i^^(u;)i.Au; + - ^ i?V(^«)(9,7fc, , 'Zj/fe, )l (3-30) 



and similarly, by Ito's formula, 

d 

2 

< \\Dipiw)\\ ■ vN^WwW + C\\D^p{w)\\ 

< CN'^exp{e\\w\\^) {\\Dip{w)\\ + ||i:'V(w)||) . 

The result for g^ is proved in a similar manner. The equality p.28p is a consequence 
of Ito's formula if G C^CHn), and a density argument proves it for the general 
case. □ 

4. Rates of convergence. We are now in the situation to prove estimates for 
the convergence of both splitting schemes and cubature methods. 
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4.1. Splitting methods. Lemma 4.1. For all k > 0, P^^ Bp {Hn) C Bp {nN) 
and supjgjQ ||P/^(^||.0-,i: < ifTllfpll ^^.A: with some constant Kt independent of (p. 

Proof. This is proved using similar estimates as those given in [13, Lemma 4.10, 
1. and 3.]. □ 

Using Lemma [3751 the method of [15] yields the following convergence estimate. 

Theorem 4.2. Let Q^^^-j '■= ^At'/2-^At'^^Ai'/2 denote the Strang splitting approx- 
imation of P^f using P^l^ and P^p ■ For any fj < rj/2, there exists C = Cts^ > 
such that for all Lp £ B'j^^ [T-Ln) and n e N, 

\\P^V - (QrT/«))>ll'A. < GrN^n-'^MU^.^. (4.1) 

Note that if 1^9 e C^(L^) is such that for some fj < rj, 

sup ?/'77(w)"^ll£'-''/f(w)||i((L2)g,i:R) < cx) for j = 0,. . .,6, (4.2) 

then ip\-H^ G bPCHn) for all N eN. Furthermore, with jj < r]/2 implies ((3T9)) . 
Thus, we obtain the following result. 

Corollary 4.3. Assume that ip satisfies ()4.2p with fj < rj/2. For any T > 

and Wo G V.'^ , there exists C = Cwo.T.tp > such that for all ri G N 

iPrfiwo) - (OfT/„))>|H„(^o)| <C{N-' + N'n-') . (4.3) 

Proof. The combination of Theorem 13.41 and Theorem 14.21 allows us to conclude 
the desired estimate. □ 

Remark 4.4. We see here an important advantage of the second order splitting in 
comparison to a possible first order splitting. There, in the second term, the instability 
would be of the order , but the convergence would only be of first order, . 
Therefore, we can choose significantly larger here while still obtaining a stable 
method. Nevertheless, we have to stress that the given error estimate is far from what 
we would expect to obtain, see also the numerical results in Section\^ 

4.2. Cubature methods. We give a short overview of cubature methods. For 
more detailed accounts, see p51 [T^ [5]. 

Fix M e N. A set of paths = (wf )^^o : [0, 1] ^ R'^+i, ujf{s) = s, i = 1, . . . , M, 
of bounded variation with uji{0) = and weights Xi > 0, i = 1,...,M, is called 
cubature formula on Wiener space of order m if 

E[[ ■■■ [ odWi' ■ ■ ■ o dWl'] 

M 

= E^' / / da.f(ti)...dc.f(t,). (4.4) 

J JO<ti<---<tk<l 

Here, {ji, . . . ,jk) £ {0, . . . , d}'' runs through all multiindices satisfying 

+ = 0} < m. (4.5) 

This means that the convex combination of iterated integrals along the cubature paths 
up to order m equals the expected value of the corresponding iterated Stratonovich in- 
tegrals along d-dimensional Brownian motion. To scale the paths to an interval [0, At], 



dw 
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we define [0,At] ^ M^+i by := t and J,'^''^^\t) V^tLof (^). 

The cubature approximations of tlie spectral Galerkin discretisation of the stochastic 
Navier-Stokes equations over a time step of size At are then given by 

jv(s, wo;a;{'^*^) = (^iyAwNis,wo;ujl^*^) + 7rArB(/Cti;iv(s, wq; wj'^*^))) ds 

+ ^,,A,dc.l^*)^^(.). (4.6) 

Here, we apply that the noise is purely additive, entailing that the Ito and Stratonovich 
integrals of the noise terms coincide. The cubature approximation of the Markov 
semigroup P^^ reads 

M 

:=J2^^fi^N{At,wo■,ioi^'^)). (4.7) 

1=1 

To prove stability of the cubature approximation, we require that the quadrature 
formula induced by the cubature scheme is symmetric, i.e., for all i = l,...,Af, 
there exists a unique i' G {1, . . . , Af} such that Xi = Xii and w- {At) = —ujI, (At) for 
j = 1, . . . , d. This induces a corresponding symmetry for oj^^*^ . Many known cubature 
formulas satisfy such a property, consider e.g. the paths given in [23 . Moreover, given 
an arbitrary cubature formula, it is easy to construct a symmetric one from it by 
adding the reflected paths. 

Our use of this assumption is to prove an estimate for the moment generating 
function of the cubature paths at At. 

Lemma 4.5. Assume that the quadrature formula induced by the cubature scheme 
is symmetric. Then, for all continuous f : M'' — > M., 

M 

j:Kfiur'^\At),...,.r^^\At)) 

i=l 

= \Y.^.{f{4''''^\M),...,u:r'Mt)) 



2 

i=l 



+ /(-u;(^*)^i(At), . . . , ~J^^'^'\At))) . (4.8) 




In particular, J2i=i ■ • ■ , uj^'^*'' (At)) — if f is odd. 

This implies 

(At)^ <e^p(^^Atpu'^ . (4.9) 

Proof. The first two claims are clear. For the estimate of the moment generating 
function, note that, as |w,*'^*^'^(At)| < CVAt and {21)] < 2^£!, 

M d oo ^ M d , 

5:A.exp(5:.,.('^*)^^(Ai)) =E^E^«(E-.-.^^'"(^^)) ' (4-10) 

i=l j=l k=0 ' 1=1 j=l 

M d ^/) ^ d 



1=0 ^ >' i=i j=i j=i 
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which proves the given estimate. □ 

Theorem 4.6. Assume that the quadrature formula induced by the cubature 
scheme is symmetric. Then, there exist % > and £ > 0, depending only on the given 
problem data, but not on the discretisation parameter N, such that with a constant 
C > independent of At and N , 

IIQ(it)/IU. <exp(CAt)||/||^^ (4.11) 
for At e iO,s], 7] e (0,r/o], and f e B^''{nN). 

Proof. Set wn{s) := ^^(s, wo; w-^*^) and V^{wn) ■= vAwn + 7rjv-B(/Cwjv,wjv)- 
For every a > 0, 

exp(as)||«;w(s)||2 - Ww^m? = r exp{ar){a\\wN{r)\\^ + 2{V'' {wr,{r)),WN{r)))dr 

Jo 

+ 2V r eMar){qjfk„WMir))dcui^'^''{r). (4.12) 
i=i 

Applying Fubini's theorem and integration by parts to 

/ exp(Q:r)dw|^*^'''(r) = exp(ar)a;-^*^'''(T) — exp(a(T)a;-^*^'''((T) 

J a 

-a [ wp^'^'(r)exp(ar)dr, (4.13) 



we obtain that 

/ exp(ar)(q'j/fc.,w;jv(r-))da;^^^*^'^(r) = (9j/fc.,wjv(0)) / exp(ar)dw|'^*''-'(r) 
Jo Jo 

+ /'exp(ar) r {qjfk,,^ {wN{q)))dqdJ,^'^'\r) 
Jo Jo 

+ J2 /'exp(ar) /''(g,/,, , <j,/fe.)dc.f *''^(g)da;f (r) 

< {Qjfk, , wn{0)) exp(as)a;|'^*^'^(s) + Cexp(as)||w;iv(0)||^Ai + Cexp(as)s2 

+ CVAfJ cxp(ag)||y^(w;v(g))||-3dg + Cexp(as)s. (4.14) 

Jo 

This yields, as {wN,V^iwN)) = -i^\\wn\\1 and ||y^(M;jv)||-3 < ||w;iv|| + C||w;^f, 
\\wN{At)\\^ < {eM-aAt)+CAt)\\wN{0)f (4.15) 

d 

+ 2 Y^iQjfk, , w;jv(0))a;f *^'^(Ai) + CAt + C{Atf 

+ y exp(a(9 - At)) ((a + C^/Ai)\\wN{q)f + C^/Ai\\wN{q)\\ - 2u\\wN{q)\\l) dq. 

Fix a = v. As ||w7v||i > ||wAr||, wc can choose e > such that for At G (0,e], 

v\\wNf + CVAi{\\wN\\ + \\wNf) - 2iy\\wN\\l < CAt. (4.16) 
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By Lemma 331 

M d d 

J2 A. exp(2r;^((?,A.^ , wn{0)}J,^'^^' (At)) < exp(rj'CAtJ2{Q,fk, , wn{0))') 
1=1 j=i j=i 

< exp(?72CAi||u;A,(0)f ). (4.17) 

Hence, for At e (0,e], 

M 

^A, eMv\\wN{At, wo; J.'^'^^'W) (4.18) 

i=l 

< exp(cAt + f]\\wN{0)f{exp{~i^At)+r]CAt)y 

Choosing ryo > small enough, we see that 

exp(-iyA<) + -qCAt < 1 for At € (0, e] and rj G (0, tjq]. (4.19) 

The claim is thus proved. □ 

Remark 4.7. It is clear from the proof that a corresponding result can also be 
shown in the space continuous case. As remarked before in the context of the splitting 
scheme, however, we are not able to derive rates of convergence in this setting, which 
is why we focus on the space discrete case. 

As it is straightforward to obtain an asymptotic expansion of by the fun- 

damental theorem of calculus (see [231 1313]), we have the following result. 

Theorem 4.8. Fix rj > .small enough. Given T > and fj < ri/2 and assuming 
that m is odd, there exist constants e > Q and C = CT,fi > such that for all 
ip e BpiHw) andneN with T/n < e, 

\\PTV~iQiT/n)rvU„ < C^iV2^(At)^||^|U„6. (4.20) 

The following result is a version of Corollarv 14.31 for cubature approximations. 

Corollary 4.9. Suppose m odd, and fix rj > small enough. Assume that ip 
satisfies ()4.2p with fj < r]/2. For any T > and wq S H^, there exists e > and 
C = Cwo,T,ip > .such that for all n G N with T/n < e, 

\Pt^{wo) - (g|^/„))>|«„(u'o)| < C (iV-i + TV^^n-^) . (4.21) 



5. Numerical examples. We consider the problem of approximating (|2.ip with 
ly - 10-2, wo = 0, d = 4, = 1, J = 1, . . . , 4, and h = (1, 0), fcz = (-1, 0), fcg = (1, 1) 
and ^4 = (—1,-1). [131 Example 2.5] shows that the dynamics generated by this 
process are ergodic. We aim to find estimates for E[||w(l, 0)|1], E[|| w(l, 0)|| _i] and 
E[||it;(l, 0)||+i]. We remark that the first and second values equal, up to a constant, 
the mean enstrophy and energy, respectively. Furthermore, control of the norm 
of w(l,0) means control of the norm of JCw{l,0), which in turn implies that we 
can take point evaluations of JCw{l, 0) due to the Sobolev embedding theorems in two 
dimensions. This is important in the evaluation of cross correlations. 

Our numerical simulations are performed using a splitting scheme, the symmet- 
rically weighted sequential splitting 

Q^^n - I + , (5.1) 
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going back at least to [211 equation (25)] and being of second order for problems that 
are smooth enough. 

We apply a Monte Carlo method. For a single realisation, we have to solve, 
alternatingly, a time-dependent Euler equation and an Ornstein-Uhlenbeck equation. 
Note that the solution of the Ornstein-Uhlenbeck equation follows a Gaussian process, 
and its distribution is therefore explicitly known. To discretise the Euler equation, 
we apply the standard RK4 scheme. While the Heun method, i.e., an RK2 scheme, 
provides the correct order such that the entire approximation is of second order, see 
[28) . it has suboptimal stability properties, leading to strong step size restrictions, 
see [51 Section D.2.5]. In this regard, see also jlTj for issues of stability of the Euler- 
Maruyama scheme for equations with non- globally Lipschitz coefficients. As we apply 
the FFT to determine the value of {ICwn ■'^)wn efficiently, we observe aliasing effects, 
which are reduced by the use of the 2/3 dealiasing, see [71 Section 3.3.2]. 

To find the expected values in the definition of P^j"^, we use quasi- Monte Carlo 
integration, applying the Sobol' sequences of Joe and Kuo [50]. Also, instead of 
simulating both terms in the definition of n' ^® ^ Bernoulli random variable 
to generate either of them, retaining the order of the approximation. 



Fig. 5.1. Error plot, increasing number of timesteps 



Figures [01 15 . 21 and 15.31 present the results of numerical calculations with increas- 
ing number of timesteps, Fourier modes, and quasi-Monte Carlo paths. All errors are 
relative, and were calculated through comparison with a reference solution found using 
K = 2^^^ quasi-Monte Carlo paths, iV = 32 and n = 128 timesteps. There, we obtained 
the values E[|lu;(l, 0)|| _i] - 1.138449630686444, E[||w(l, 0)|1] ~ 1.319968848291092, 
and E[||w(l,0)||+i] - 1.620419847035606. In Figure [O we chose the other param- 
eters to be if = 2^^ and N = 32; in Figure |Ll K = 2^^ and n = 128; and in 
Figure |01 TV = 32 and n = 64. 

We clearly see that mainly the number of quasi-Monte Carlo paths limits the 
attainable accuracy. Nevertheless, with 2^^ — 4096 paths, we obtain a relative error 
of less than lO^'^, and that calculation took approximately 60 seconds running on 16 
cores of a Primergy RX200 S6 spotting 4 Intel Xeon CPU X5650 processor, each of 
which provides 6 cores. In Figure [STTl we observe that we obtain a rate of convergence 
of about 2.5 for the norm with respect to the number of time steps, which is even 
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Fig. 5.2. Error plot, increasing number of Fourier modes 
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Fig. -5.3. Error plot, increasing number of quasi-Monte Carlo paths 



more than the theoretically predicted rate of 2 and seems to result from the fact that 
we compare with numerical estimates instead of the exact value. The solution of the 
model problem is smooth (see also US] in this regard), and indeed, Figure exhibits 
spectral convergence in the number of Fourier modes. 

6. Conclusion. We have introduced and analysed novel high order approxima- 
tion schemes for the stochastic Navier-Stokes equations on the 2D torus. We prove 
high order accuracy in time and give precise estimates for the dependence on the order 
of the spectral Galerkin discretisation. Using high order cubature paths, it is possible 
to attain convergence of arbitrary order in time. 

From a practical point of view, the splitting schemes presented in this work have 
the important advantage that well-tested and robust solvers for the deterministic 
Navier-Stokes and Euler equations can be reused. Furthermore, the algorithm makes 
increasing the dimension of the driving Brownian motion easy. Numerical examples 
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establish the apphcabihty of the method to some simple, but relevant functionals. 
Appendix. Proof of Proposition [273l 

Lemma A.l. For N - Af{{), 1), j = I, . . . ,d, and S , A, B eM. with C € R small 
enough, 

E[exp(C(52+25^BA + {BNf))] 



(l-2CB2)i/2 1-2CB^J 

Proof. A direct calculation yields 
IE[exp(C7(S'2 + 2SABN + (BNf))] 

^ J^eMCiS' + 2SABy + exp (^~\y^^ dy 

1 [ ( I, ^ 2CSAB V\ , 

1 // \ „ 

which proves the result. □ 

Corollary A. 2. For independent Nj ^ A/'(0, 1), j — l,...,d, and S, Aj, 
Bj e R with C e M small enough, 

d d 

E[exp(C(52+ ^ 2SA,B,N, + Y.{B,N^f))] 



1 (ii + s-i 



2CAjSj 



2CB| 



CS-^ . (A.2) 



Proof of Proposition \2.3\ Note that 

w^{t,wo)^exp{tei^A)wo+ exp{{t - s)eiyA)^qjfk^dW^ . (A.3) 

Denoting by Aj the eigenvalue of fk^ with respect to the operator evA, evAf}^. = 
Xjfkj , we see that 

/ exp((t - s)eiyA)QdWs = / e^P((* " «)-^j)9j7fe.dM^i . (A.4) 
Jo Jo 

The coefficient := Jjj*exp((i — s)Xkj)dW^ is normally distributed, more precisely, 
Z{ ~ A/"| 0, ^~°^P(2*^'°j) y jn particular, with ^(i) := exp(te:yA), 



e^v(^r^\\S{t)w + J2q,Z{fkX 
i=i 



(A.5) 
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Note 

d 

\\s{t)w+Y.'i=zih^r = \\s{t) 



3 = i 
d 



and apply Corollary |X2] with C = 77, 5 = \\S{t)w\\, Aj = \\s{t)wty\\q^h ■ llo ^^'^ 

\ 1/2 

l-cxp(2tA,^.)\ . , .2 



- Il9.7fcj| ii;-^ j . As < 1 and 

1 - 2CB| = 1 - 2r,\\q,h^f^—^^:^B^ > expi2ujt) (A.7) 

— 2Afe^- 

for > a; > Afc^ and < 77 < and, similarly, 

for a > and n < min,-=i d tj 1 ^ 1?" j — nr: we obtain 

P^ipniw) < exp{—dtuj) exp{rj\\w\\'^) — exp{~dtuj)^lj,^{w), (A. 9) 

the required result. □ 
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